// PANERI Jérémy
// PEREZ Julien
// M1S8 Info
// Projet CUDA

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#include <cuda.h>
#include <cuda_runtime.h>
#include <helper_cuda.h>

const int nb_p = 64*1;
const float dt=0.05;
const float G=6.672e-11;


// Fonction CPU
void calcul_cpu(float4 *X, float3 *V, float3 *A) {
	int i;
	//calcul position
	for(i=0;i<nb_p;i++) {
		X[i].x+=V[i].x*dt;
		X[i].y+=V[i].y*dt;
		X[i].z+=V[i].z*dt;
	}
	
	//calcul vitesse
	for(i=0;i<nb_p;i++) {
		V[i].x+=A[i].x*dt;
		V[i].y+=A[i].y*dt;
		V[i].z+=A[i].z*dt;
	}
	
	//calcul acceleration
	for(i=0;i<nb_p;i++) {
		float3 a={0,0,0};
		float4 x=X[i];
		for(int j=0;j<nb_p;j++) {
			if(i!=j) {
				float r=sqrtf((X[j].x-x.x)*(X[j].x-x.x)+(X[j].y-x.y)*(X[j].y-x.y)+(X[j].z-x.z)*(X[j].z-x.z));
				r=r*r*r;
				float s=X[j].w/r;
				a.x+=s*(X[j].x-x.x);
				a.y+=s*(X[j].y-x.y);
				a.z+=s*(X[j].z-x.z);
			}
		}
		A[i].x=G*a.x;
		A[i].y=G*a.y;
		A[i].z=G*a.z;
	}
}


// Fonctions d'affichage
void afficheX(int nb_p, float4 *X, int t)
{
	int i;
	printf("Itération %d\n",t);
	for (i=0;i<nb_p;i++) {
		printf ("%d : %f\t%f\t%f\t%f\n",i,X[i].x,X[i].y,X[i].z,X[i].w);
	}	
	printf("\n");
}

void afficheAV(int nb_p, float3 *AV, int t)
{
	int i;
	printf("Itération %d\n",t);
	for (i=0;i<nb_p;i++) {
		printf ("%d : %.15f\t%.15f\t%.15f\n",i,AV[i].x,AV[i].y,AV[i].z);
	}	
	printf("\n");
}




// Calcule la position sur GPU
__global__ void calculPosition (float4* X, float3* V){
	int i = threadIdx.x + blockDim.x*blockIdx.x;
	X[i].x += V[i].x * dt;
	X[i].y += V[i].y * dt;
	X[i].z += V[i].z * dt;
}

// Calcule la vitesse sur GPU
__global__ void calculVitesse (float3* V, float3* A){
	int i = threadIdx.x + blockDim.x*blockIdx.x;
	V[i].x += A[i].x * dt;
	V[i].y += A[i].y * dt;
	V[i].z += A[i].z * dt;
}

// Calcule l'accélération sur GPU
__global__ void calculAcceleration(float4* X, float3* A){
	int i = threadIdx.x + blockDim.x*blockIdx.x;
	float3 a={0,0,0};
	float4 x=X[i];
	
	for(int j=0;j<nb_p;j++) {
		if(i!=j) {
			float r=sqrtf((X[j].x-x.x)*(X[j].x-x.x)+(X[j].y-x.y)*(X[j].y-x.y)+(X[j].z-x.z)*(X[j].z-x.z));
			r=r*r*r;
			float s=X[j].w/r;
			a.x+=s*(X[j].x-x.x);
			a.y+=s*(X[j].y-x.y);
			a.z+=s*(X[j].z-x.z);
		}
	}
	
	A[i].x=G*a.x;
	A[i].y=G*a.y;
	A[i].z=G*a.z;
}


// Fonction principale d'appel GPU
void calcul_gpu(float4 *X, float3 *V, float3 *A) {
	calculPosition<<<nb_p/16,16>>>(X, V);
	calculVitesse<<<nb_p/16,16>>>(V, A);
	calculAcceleration<<<nb_p/16,16>>>(X, A);
}




// Programme principal
int main( int argc, char** argv)
{
	int t = 0;
	const int t_max=50;
	float time;
	
	clock_t startCPU, stopCPU;
	cudaEvent_t startCUDA, stopCUDA;


	// Création et initialisation des variables CPU
	float4 X[nb_p];
	float3 A[nb_p];
	float3 V[nb_p];
	
	float4 Xbis[nb_p];
	float3 Abis[nb_p];
	float3 Vbis[nb_p];
	
	for (int i=0;i<nb_p;i++) {
		X[i].x= 1000*(rand() / (float)RAND_MAX);
		X[i].y= 1000*(rand() / (float)RAND_MAX);
		X[i].z= 1000*(rand() / (float)RAND_MAX);
		X[i].w= 100*(rand() / (float)RAND_MAX);
		A[i].x= 2*(rand() / (float)RAND_MAX);
		A[i].y= 2*(rand() / (float)RAND_MAX);
		A[i].z= 2*(rand() / (float)RAND_MAX);
		V[i].x= 5*(rand() / (float)RAND_MAX);
		V[i].y= 5*(rand() / (float)RAND_MAX);
		V[i].z= 5*(rand() / (float)RAND_MAX);
	}

	
	// Création des variables CUDA
	float4 *cudaX = NULL;
	float3 *cudaA = NULL;
	float3 *cudaV = NULL;

	// Allocation mémoire CUDA
	cudaMalloc((void**)&cudaX, sizeof(float4)*nb_p);
	cudaMalloc((void**)&cudaA, sizeof(float3)*nb_p);
	cudaMalloc((void**)&cudaV, sizeof(float3)*nb_p);

	// Copie du CPU vers GPU
	cudaMemcpy(cudaX, X, sizeof(float4)*nb_p, cudaMemcpyHostToDevice);
	cudaMemcpy(cudaA, A, sizeof(float3)*nb_p, cudaMemcpyHostToDevice);
	cudaMemcpy(cudaV, V, sizeof(float3)*nb_p, cudaMemcpyHostToDevice);


	// Exécution CPU
	startCPU = clock();
	for(t=0;t<t_max;t++) {
		calcul_cpu(X, V, A);
	}
	stopCPU = clock();
	
	printf("\n- Temps CPU : %fsec -\n\n", (float)(stopCPU-startCPU)/CLOCKS_PER_SEC);
	afficheX(10,X,t);
	afficheAV(10,A,t);
	afficheAV(10,V,t);


	// Exécution GPU
	cudaEventCreate(&startCUDA);
	cudaEventCreate(&stopCUDA);
	cudaEventRecord(startCUDA, 0);

	for(t=0;t<t_max;t++) {
		calcul_gpu(cudaX, cudaV, cudaA);
	}
	
	cudaEventRecord(stopCUDA, 0);
	cudaEventSynchronize(stopCUDA);
	cudaEventElapsedTime(&time, startCUDA, stopCUDA);

	// Copie du GPU vers CPU
	cudaMemcpy(Xbis, cudaX, sizeof(float4)*nb_p, cudaMemcpyDeviceToHost);
	cudaMemcpy(Abis, cudaA, sizeof(float3)*nb_p, cudaMemcpyDeviceToHost);
	cudaMemcpy(Vbis, cudaV, sizeof(float3)*nb_p, cudaMemcpyDeviceToHost);

	printf("\n- Temps GPU : %fsec -\n\n", time/1000);
	afficheX(10,Xbis,t);
	afficheAV(10,Abis,t);
	afficheAV(10,Vbis,t);


	// Libération mémoire
	cudaFree(cudaX);
	cudaFree(cudaA);
	cudaFree(cudaV);

	return EXIT_SUCCESS;
}

